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SUMMARY 


This report describes a program for the numerical solution 
by an explicit finite difference technique of the momentum, 
species conservation, and energy equations for an equilibrium or 
chemically reacting air boundary layer. Laminar and several 
turbulent viscosity models can be applied through the thickness 
of the bovndary layer. Seven air species are included: O25, O, 
Ne, N, NO, NOt and em. 


Laminar and turbulent effects in any of the equations can 
be selectively controlled by input indicators. The various 
viscosity models can be applied to the entire thickness of the 
boundary layer or the boundary layer can be divided into two 
segments and different viscosity models applied to each segment. 
Provision has also been mace for boundary layer swallowing of 
new inviscid streamlines. Properties along the inviscid stream- 
line are computed to establish the varying streamwise outer edge 
conditions. 
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aE INTRODUCTION 
This program is written in FORTRAN II language for the 
IBM 7090-4. [It solves the non-equilibrium and equilibrium air 


boundary layer conservation equations for momentum, energy and 


species mass fractions. These parabolic partial differential 
equations are solved in the domain defined by the body surface 
and the inviscid edge nf the boundary layer. The technique of 
solution is an explicit finite differerice method. The boundary 
conditions at the body are specified as input data but those at 
the inviscid outer edge are not known a priori and are therefore 
computed as the calculation progresses. Edge conditions are 


determined from the solution of ordinary differential equations 


which yieia the velocity and species along an approximate 
inviscid streamline from the shock to the body station of 
interest. Initial conditions to start the problem require the 
distributions of the dependent variables of velocity, total 
enthalpy, and species mass fractions. 

A number of significant capabilities are provided by the 
program. The multicomponent air species considered are Os, O, 
Ng, N, NO, No* and e. Laminar and several turbulent viscosity 


representations may be optionally selected for the diffusional 


transport terms applied to any of the partial differential 
equations. The boundary layer thickness can be divided into 
two segments and different viscosity models applied to each of 
the segments. Various two-dimensional or axisymmetric body 
geometrical configurations are possible. Similarly, shock shapes 
consisting of conical and parabolic segments are acceptable. 
Input tolerances upon the slope at the outer edge of the enthalpy, 
velocity or species profiles govern the boundary layer swallowing 
which results from the growth of the viscous layer into the 
inviscid flow. 

The program described herein is based upon the analysis 
reported in Part I. The engineering definitions and the prepara- 


tion of input data for application of the program to given boundary 


layer problems are discussed in the Part III manual. 

Typical hypersonic reentry body configurations have been 
investigated. Computer time required for their solution ranges 
from 3/4 to 4 hours, depending on body length, altitude, free 


stream velocities and mode of computation. 
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II. BOUNDARY AYER EQUATIONS 


In this section the explicit finite difference method solu- 
tion of the boundary layer conservation equations is described. 
A. Coordinate Systems 
The program continuously calculates in four different 
coordinate systems. Thus for the permissible body and shock 
geometry configurations, the following systems shown in the 
sketch below are defined: 

l. Overall - a basic reference system designated as 
C1, Cs. The units of these coordinates are feet. 

2. Physical - an x-y system where x, in feet, is along 
the body and y, in feet, is normal to the body. 

3. Body and von Mises stream function - an x-W system 
in which the finite difference solutions will be 
solved, where x, in feet, is along the body and ¥ 
is the transformed normal coordinate and represents 
the stream function defined by W = [purfay. For 
€ = 1, Wis in slugs/sec; for ¢ = a ieee ry 
slugs/ft-secs. The sketch below shows the grid in 


the x-V¥ plane with the paths of the numerical step 


indicated. 


Inviscid streamline - a fourth system is a one- 


dimensional coordinate, s, along the (approximate 


inviscid) streamline. The origin of this coor- 


dinate system is the shock streamtube intersection 


point. 


Ce 


Ci 


Coordinate Systems for a Blunt Body 


m=~ 


Finite Difference Grid 


B. Momentum Equation 
1. Analytic form 
’ A 
ou ee OD ee ek nol 
ox pu ax ay 
Seo es 
T = pb Sy pur 
2. Finite difference form 
€ 
u = - ao lel r ' a fh | 
m+l1,n m,n oe . Un, es 6 m | m,n+1/2 m,n-1/2) 
Note: Near the wall a series solution has been 
assumed since u~ 0. The expressions used are; 
= Os 
u=a,V + b, ¥ 
i 5 
c. =c, ta. Wtb, ¥ 
1 al 1 in 
Ww 
Ys 
H=H + agW° + bs ¥ 
Ww 
L 
pu=alv?+b,v 
Using the variables evaluated at AW and 2AY, the 
ais and b's are determined. Thus from the 
‘ AW 
series solution the values at = required in the 
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€ 
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4. 


finite difference solution of the momentun, 
species and energy equations are then 


evaluated. 


The difference in shear stress 


€ L( L L 
mee ee dn n a Dene targl Ym, n+1 “nn 
2 2 1B, n+1/2 AV 
€ Tifa L eee 
Pao mene el eae | gato? “m,n -m,n-1 
2 2 *"m,n-1/2 AW 


Viscosity models 


5 


3.04566 x 107° 9" 
se ents ey = lems ae 


’ 


aed kl Tet? Y big T 74 
Bio + Bah 6+ Bs 6+ Bab 


b. ‘Turbulent ie 


t Bey? + Bout” + ey ra where the B's are 


input coefficients and the Age are defined as 


foiliows: 


i. Law of the wall - dimensional coordinate, 
2 
p +p 
ee aly + AY 375 aL25 (Fx)? { m,nt+tl m,n 
m,n+1/2 m,n aE) (ous) 2 
( e m,n m,n+ 
) 
(pur *) +(pur *) u -u 
m,nt+l m,n ntl m,n 
2 AY 
2 
p +p 
yit aly AW wee E AWE (FK)? mn m,ntl 
m,n-1/2 m,n (pae®) ( ar’) 2 
m,n=-1 ie m,n 
€ € 
ous Im nee ) mp n=1 m,n -m,n-1 
2 AY 
ii. Law of the wake - dimensional coordinate, 
T2 - (rmentl mn} {-m.ntl" nn} (yu) 
E Wa, nt1/2 iS 2 2 J ee'm 
+ 
T2 ES (aaa (hun nna oy oa) 
Wo, n-1/2 2 2 e e’m 


where K(N) is an input function of V 


K(N) =a+b¥+c¥ 
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Ye° 


iii. Laminar contribution 
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iv. Eddy viscosity fit - dimensional coordinate, y 
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v. Law of the wall - streamline coordinate, YW 
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vi. Law of the wake - streamline coordinate, , 


ks _ fess 7au esa ue 
m,n+1/2 — 2 2 € 
(Por). 
ae oa Fanner Facey ue 
m,n-1/2 — 2 2 € 
(Por). 


vii. Eddy viscosity fit - streamline coordinate, 


W V+ & vt Ae 
-14- .08 |{1- 7 - .06 |{1- TY = 
(Puor ) e e 
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Boundary layer swallowing and 
expansion of the boundary layer 


An input number qa is used as a test value for 
calling a new inviscid streamline and to increment 
the ¥W grid. Two options are available based upon the 
value of eo If = 0, the streamline test is not 
made. If ce is set as any positive number (used as 


the tolerance), the test 


ee a eee 


S € 
u 
Unty (htL) 


is made. If the test fails, a new streamline is 
computed, or constant edge conditions are carried 


over. When the incremented expansion of the W grid 


u_ (bd) 


ube 1) eee 


is called for then set Wey (b-9) = 


= = = ‘ 
for J 2, 3, 4, and Wey 1) Ut (ht2;- 


Boundary conditions 

a. Initial distributicns - an input velocity profile 
vs. physical coordinates. 

b. Inner value at the wall - u = 0. 

c. Outer edge - velocity at the outer edge is 
obtained from the finite difference equation 
of “B" or the inviscid streamline solution of 

du 


U = — 


ae 22 accordingly as the momentum, species 


Dlr 


or enthalpy slope tests are satisfied or not 


(see Subsections II.B.5, I3.c.5, II.D.5 for 


formulation of tests). 


d. Inviscid streamlines - two inviscid streamlines 
external to the boundary layer which represent 


the locally uniform (assuned) external flow are 


always computec. 


C. Species Equation 
1. Analytic form 
EAP EU ie ee oan 
ax pu aw’ nae 
2. Finite difference form 
A x(w. A 
= (c.) Nie are , AX oe | = 
mtl,n “Enea (Pu) AW ~m,n m,nt+1/2 m,n-1/2 


See Section II.B.2 for note reg :ding solution near 


the wall. 
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3. The species diffusion terms are given by: 


€ € LL LL 
WA Wes eB Megas yn n L H Dig [u PEA 
= 1/2 2 ToL Z | ae 
Ss _D SD 
Cc m,ntl Cc m,n 
i 
T T 
cet he D3 | pa aleeek oy nen CAU yaa 
SH nt1/2 T 7 T T AW 
S D S D 
c m,nt+l 2 m,n 
i 
€ € L L 
+ é : : 
(Pur )mapn fae mn L uw Di uw met 
Se 2 po a eer: Deo 
S D SD 
e m,n C m,n-L 
a 
watt = Dei : (2 Dei (Cia nT (oa nal 
3 — 
m,n-J./2 a a st pt AY 
Cc m,n Cc m,n-l 
4. Chemistry production terms 
a. Formulation 
w, (W,x) = w.(T,P,C,), 5 aes Senay 


Oe 


Se 


Se 


N2 


where 


= 
M 
=e -A+t+— p 
Myo 
M M 
p 
= A + =° (c-p-z) - = F 
NO Not 
p M 
Myo 
M 
p 
= I: ete (C+D+E) - A 
L NO not j 
p 
== [-c-pie-x] 
bt 
=F [rox | 
M 
eo [Fes] 
Not 
eoupye C eee 2 | 96573302 x 1028 
ra P ao02® see ou lacs sit 
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Cc . -D.__/KT 
ae 2L at 1 N2 as 2 26 
3 : ma Ki p P anc o® PCy 9.6573302x 10 
N 1 
ante C fi ~Dyo/KT 8 
bat Ree ees = 2 
Cc MM Y v7 KP Paxoo® PCA 9.6573302 x 10 
NO 1 
. Pno7Po2 
_ —_NG_ KT 23 
D 1M Kg? for CNC o® CO2 CN 3.1076201 x10 
N ©2 
P27 no 
EK Ae KT 
ee Oona - - 23 
E M. PY Cy2So® CN0CN 3.1076201 x10 
Tye No 
LK 
: Fe eke p K C.Cc.e <r - (Cc. +)c 3.1076201x 10? 
M. O-N NO e : 
€ -I_ /KT 
L b eee 2 NO 26 
= aa - + ; 0 
K M. : M, KK p P ae No® P(CLo Ce, 9.6573302-% 1 
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ET 2274 = 1 - @ 2274/7 
ER o7aGeest sets 
Erogsos = een 
ET 1130 = 3 + 2 e 11300/T 


Then 
M Ye 2 a 
¢. = —s (nm KP) (4-152 (ev 2274) (4E22 228)" | 3 97871 x 1072? 
ac 27h. Oo Ta} ET 1130 
8M : 
2 : a 
o = —S (nm KT) (3 789! (er 3395) 3.87871x 107? 
aN h T 
p 
a My S/o 
fe) fa 4.906 ETT 228 =) 
o eh coaaece — Bie eS a 
daNo nip (2m) KT) be [4298 42308) (en 2740) | ET 178 3.87871 x10 
8/5 
o= (faa ee ET A ET 1130 | 
MoM 4.161] \ET 2274] | ETT 228 ET 178 
S/o 
esate Myon 5.789\ /pr 3395 ET 178 
MoM 4.906] |EeT 2740] l|ETT 228 
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K Myot) “Me T 1 | 
mm 
MM, 5.789 Br 3395 ETT 2298 
M (2 xn)?” 
™ 
meen =: e ET 2740 2 = 
P se Secu eee (.4237) [BE 2240 ET 178 3.87871x 10 
b. Input Az, = 1 for reacting system 


= 0 for frozen system; bypass all ws 


= 1 for equilibrium system. 


MW 


10 for reacting system with equi- 
librium at the wall 


5. Equilibrium mass action laws 
As indicated in II.c.4.b, an input of -1 for 
Az, activates the logic for obtaining the species 


from the following equilibrium relationships: 


B/> ple iets = ei = = = 
908.19 T 2 0S . BE OO Se eheecalors 


+ 
enh ce Os 227 
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_ 228 Dee 29501 
Ye T T Gata tee. -9 
e (0) = 2237-7 1.6 e + .2 e e A(1.386x10_) 
m+1 p 
5/2 x 
© (Ne) = 228-39. B® (2.428 x10 a) 
mie ey ees | 
Pp \l-e 2 
; 27698 41520 56544 
Ye ei = as = 
PAGOT! Ges) nen, en ce : a, 02 He 6 2" Satria don 
oO — ————  — oe 
F | _ 178) _ 11063 
a mee -16 
ss 5) = 23414 tT Ite plietews BEAD (On GoclOles) 
m+1 (NO) = 2740 
p ke A: 
epevalels 
-11 Ly 
veer .1597x10° Te Cree) ce tN) 
“m+1 Sry _ 3395 _ 228 a 
l-e z 1.6 e f + .2e ¥ 
; 
(NOt) 


(e-) = ——>—- c 
m+l (233.9)? m+1 
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Ac 


>| 


BT 


AC 


A and B are obtained from an iteration pro- 
cedure which is satisfied when two successive 


values of B agree to within 


1 
ee ©} 
1000 C 
The iteration procedure is initiated with B= 


and then the following calculations take place in 


the sequence shown: 


a) | jigeyse 
3/2 ie = 8/s fm ~- 
D\e2o7 teh, were ee AO MOsan Aen nyecumen el aT Ons 
Es 
B/, =o 
_ (908.19 T ”)B,P(1338.7x10~) 
2 = 
(eee (BT) *+Ac hy 
Sh DRE tee] VAeraraee Rene 
(908.19 TB, ) (4x10?) 
56544 11063 
Ya a, T - 6/2 - T -18 
_ \1461.9 T Bs e A\ 934.147 Ey e 10 
= = 
Bf “le 
(178.39 T )p(5060.8x10  ) 


Eg 


ee / (BT)? +Ac les 
a a -12 
Ci Oe39.T yao 10m |) 
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Boundary layer swallowing and expansion 
of the boundary layer 


An input number 2 ls used as a test value 
for calling a new inviscid streamline and to incre- 
ment the W grid. Two options =:e available based 
upon the value of eae i ie EQ = 0, the streamline 
test is made. If Eu is set as any positive number 


(used as the tolerance) the test 


DS ps ait, 


is made. If the test fails, a new streamline is 
computed, or constant edge conditions are carried 
over. When the incremented expansion of the V 
grid is called for then set Cnty (EI) 


c(h-d) 


= ¢_ (il) Cty (ht2). for J = 2, 3, 4 and 


cy iol) = Cty (Et) - 


Boundary conditions 
a. Initial - C; (y) profiles are input for all spe- 


cies except Np and e, which are computed from: 


Mee : 
SiN) =e a (NOT, 
sen ) Moot ( ) 
cy (N) = 1- Be, i # No 


2 


XxX x 
b. Wall - Cty (2) = C, (I) + Cs (I) Rn + e3 (1) [2 


3 


4 
26 e4 (2)(2] + Cg (I) 


Rn where the C, through 


Cs are input for each specie except e. 


fualisiuas 


G 


Outer edge - Species at the outer edge are 
obtained from the finite difference equation in 
Subsection II.C.2. or the inviscid solution of 
dc, w. 

2 : 
Mo as py accordingly as the momentum, species 
or enthalpy slope tests are satisfied or not 
(see “Subsections 1108.5, 11.C.5., £2.05. for 
formulation of tests). 

ad. Along the two external inviscid streamlines the 

species are obtained from the inviscid equation. 

D. Energy Equation 


1. Analytic form 


2. Finite difference foim 
e Ax H H 
elem a en ~ “m AW (eS ecu, Se aaa! 


€e Ax C C 
ae ety = eye 


See Section II.B.2. for not? regarding solution near 
the wall. 
‘ Sree (ou 
3. The difference of enthalpy-kinelic energy terms 


is given by: 
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4. The energy transport due to species gradients C° 


is given by: 


ome ten Re a Beet aly ae 
4 We iL Alacer oh : oh 
e/m,nt+l e/m,n cflm,nt+l c/m,n 
Diels L 1 
(Pur) ney Ur) | [AS (Ha nad Leo) rile 2 (h +h ) 


L L (c.) -(c.) 
me Dei ‘smatl 7 mol 2} f, 2 ete 
pu m,n+l pu m,n Ay : L Ta 
e/m,nt+l e/m,n 
as lL. € Oe T_T 
a . T (eur) aoe se n Be nei 
Ss 
c/m,ntl c/m,n 
L De aa ous he a n 
a u ith +h ) - + aa Sai ~PARTC2 
= m,nt+l m,n D m,n+l \D m,n 


5. Boundary layer swallowing and expansion 
of the boundary layer 


An input number en is used as a test value for 
calling a new inviscid streamline and to increment 


the W grid. Two options are available based upon 


the value of ea If Gy = 0, the streamline test is 
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not made. If ey is set to any positive number 


(the tolerance) the test 


- hip) 
Hy (ht) Hu 1} 


H yy (ht) H 


is made. If the test fails, a new streamline is 
computed, or constant edge conditions are carried 
over. 

6. Boundary conditions 
a. Initial - an input enthalpy profile vs. physical 


coordinates. 
2 


b. Wall - Ha Rn 


+1 


= HH, + HH, || + HHS 


3 4 


Hid where the HH, through 


+ wy (2 


Rn 
HHs are input. 

c. Outer edge 

Enthalpy at the oute~ edge is obtained from 

the finite difference equation in Subsection 
II.D.2. or set equal to the input free stream 
value accordingly as the momentum, species or 
enthalpy slope tests are satisfied or not (see 
SuUbSePCEIONS Lb iabeoapy bikes, LL eve afore 


formulation of tests). 


Lit. 


INVISCID FLOW EQUATIONS 


A. Governing Equations for One-Dimensional Streamlines 


1. Momentum 


inv iste) 
u FE =_l 
inv ds P. ds 
inv 
2. Species 
dc, 
inv Ww 
ou = 
inv ds p. 
inv 


3. Total enthalpy (energy) 


Hae} 
il 


constant 


B. Tracing the Streamline from the Shock 


Reference: NACA Report #1135 - “Equations, Tables and 
Charts for Compressible Flow 1953" 


1. Intersection of streamline and shock, see sketch. 


1, 
a | omen 


2 pi uit 
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DL 


z € aw € 
J € |(27) y cos 0, ie + (2m) V, 
fe) 


90 - a, for arc, ogive 
QO, cylinder, flat plate 
cone angle, cone 


angle which the current point on the body and 
the coordinates of the center of the arc make 
with (, axis, measured in a clockwise direction 
from (C,. 


C2 


C3 


Coordinate Systems for a Blunt Body 


; 2: 


a. 


De 


Petes 


| 
| 
| 
| 


Conditions downstream of the shock 


Following will be input 


Pe 
stag 
h h 
foe} [oe] 
needed for T =— 
) c 
c 
loam Po 
De 
altitude 


y (specific heat ratio) 


- used for checking against calculated 
u in the one-dimensional streamline 
equation in order to keep u > 0. 


“TOLIM 


Following will be program set from input 


u, = u(L) 
oie = c(L,I1) 
H = H(L) 


The program calculates: 


L 
> 2 
[4(M,* sin? @-1)(y M,”* sin® 6+1)] 
ds [(y+1)* m,* sin’ @] 


a2 


where 
u 
eo 
M, = =~ Mach # 
a 
fo°) 
a, = speed of sound = 49.9 /1.8 ae 
hy 
ig eh © 
ea G K 
Pp 


6 = $~ shock wave makes with C, 


foe) 


: E M,* sin® 9 - I 


veer a ye 


dad. The program outputs: 


(a) Pos 
(b) Yds 
(c) Cia, 


3. Completion of streamline calculation 


The program sets: 


‘ = = 
U1 Sh) Une] (Ett) Ut (Lt2) 


Cty het) = C ey Ett) = C ey Et2- 2) 


IV. INPUTS 


Under this section the various input parameters are dis- 
cussed from a programming viewpoint; see Manual, Part III, for 
input data preparation. In Subsection A through C, certain 
critical parameters are discussed, and in Subsection D the input 
Gefinitions and card formats are specified. Also included in 
the latter sections of this report is a sample problem which 
illustrates the nature of the inputs. 

A. Body Geometry 

l. Classification 
Bodies will be classified accordingly as the first 
geometrical region is: 
a. Blunted As #4 0 
(a) arc (hemispherical nose) 
b. Pointed As = 0 j 
(a) cone 
(b) ogive 
(c) flat plate 
where As = standoff distance from shock to stagnation point of 


body along centerline (see IV.B.3). 


UMEEETTCESEPTSRTPTOE ALY TSE OPE TORTI 
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sections are 


Subsequent geometric regions of the body 


The following geometrical descriptions for body 


available: 
a. arc 
bh. cone 


2 dimensional or axisymmetric 

c. cylinder depending uponinput €, (or EPSI in 
the input list, Subsection IV.D) 

d. ogive 


e. flat plate 


Calculation of radius of the body and the (,, ls. 
coordinates 


a. For each region the program reads in 
(a) INDP = pressure type l = cosine 
indicator 2 = polynomial 
(see p.39, Sectionc) {| 3 = quotient 
(b) INDR = geometry indicator l = arc 
2 = cone 
3 = cylinder 
4 = ogive 
5 = flat plate 
(c) INDLR = last subregion of l = yes 


last region indicator {0 = no 


(d) INDLSR = last subregion of any 1 = yes 


region indicator 2 no 


(e) OGIVE H 
coordinates of center of arc 
with respect to C,, Ce origin 

OGIVE K 


(f) p(J), J = 1,8, pressure coefficients 
(g) cc, (K), K = 1,5, Cre: coefficients 


(h) Rn = nose or reference base radius 


cumujative upper ~ limit of this 


(i) xi = 
subregion 


(j) CONEAN = cone angle/2 (vertex semi-angle) 


B. Shock Geometry 


1. Description of the Shock 


The shock consists of up to 10 sections of the form: 


TERETE as rE eR ee TTT TM Moin: Mimi ini ii i iin 


a. Ca = aoa a be straight line 


~ = 2 
ioe OQ - an = b (Ce seo) parabola 


Type 6 is a parabola whose coefficients an! bo and c,, are specified 
as input. Type 7 is a parabola whose beginning and ending slopes 
and the €: projection of the end point of the geometrical shock 


reyion are specified as input. 


The present version of the program permits para- 
bolas to follow each other but they may not follow straight lines. 
2. Shock input is specified by giving 
a. Curve type 


b. Slope and cumulative €,; limit for a straight 
line 


c. Coefficients and cumulative ©, limit for a 
parabola or 


d. Beginning and ending slopes and cumulative 
Cy limit 


3. Section 1 of the shock has optional provision (type 7) 


for internally computing and establishiug the shock 


a es 2 
Cx ae b (te c.) 
with 
=c =0, pee eens 
n n n 14+ 48 Rn 
Rn 
where As= a Rn (physical standoff distance). == is the 
stag ay 
p 


dimensionless standoff distance. 
4, The program has provision for specifying several shock 
shapes, e.g. cOnical, parabolic and normal. Moreover optional 


modes for specifying a particular shape are also possible, e.g. 
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parabolas can be specified by inputting the coefficients of the 

analytic expression or the beginning and ending slopes. The 

following indicates the formulae for obtaining the € coordinates 

and the balances of the shock soefficients fer all permissible 2 


combinations of shock shapes. 


a 

: a. Region 1 is a cone 
3 be 
: j 
Ine cay oy = “ae inoue | 

e : : 

oie Coe ee Ce 


b. Region N-l is a parabola 


Region M is a parabola - Type 6 


Hee vailabl 
i a N-1! Cue are available 


c.. are input 
N 
se ee 


©) (1) Teenie oo oNed (le (1) "Cnea) 


c. Region N-l is a parabola 
Regicn N is a parabola - Tyre 7 
ie 


aye1! Dew Cyd are available 


C 


are input 


IN 
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7 is DI 
N a = me 
2 
; : 2m (a ne ayy) 
See SAN RNG 
N (N-1) M4 
1 


<S a 
C2 (ya) = Aner + Py-a 6S? ea) 7 S-2)? 


ad. Region N-l is a parabola 


Region N is a cone 


1 anye1' Pyed! Cured are available 
Mo is input 
ii. ¢ See 
e 2 os be oe N-1 
(N-1) 2b ay N 


= eS 2 
G2 (N-1) = eye ona [be (y-1) Cua! 


By = ©2 ya) 7 Ay St (yen) 
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e. Region N-l is a cone 
Region N is a cone 
i by are available 


Aye1' Pyeas § (wen) 


mis in,ut 
N t 


de Ca ia) = Ayer Seay * Peed) 


By = Ce cgay 7 AyySt y-2) 


Cc. Pressure Distrikution 
Associated with each body subregion is a distribution 


which is one of 3 possible types: 


ey oe=prcos: = (Newtonian) 


2 


= + ee 
Pa P3 Rn Pa Rn 


tu 
ue) 
| 


Py 


where Py through p, are inputs. 


40 


SYMBOLS 


*LP2 
NS 


NPSI 


INDSTR 
FNDSSL 
FNDSSH 


IALT 


c(N,I) 


*RESTAR 


Input Definitions, Card Formats and Symbols 


DEFINITION 
# Wpts + 2 
# species 


print interval density in W direction e.g. 
every point, every second point, etc. 


streamline starting point indicator 

ower limit for step size control 

upper limit for step size control 
altitude - in kft 

laminar indicator for momentum equation 
turbulent indicator for momentum equation 


chemistry indicator for equilibrium or non- 
equilibrium calculations (see p.18, Section b) 


laminar indicator for species equation 
turbulent indicator for species equation 
laminar indicator for energy equation 
turbulent indicator for energy equation 
laminar indicator for energy equation 
turbulent indicator for energy equation 


species mass fraction profiles N=1,2...LP2; 
iL =1,2,4,5;6 


restart indicator; l = original run, 
Q = continuation run 


FORMAT 


E 10.3 


E 10.3 


EtOnS 
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SYMBOL TEFINITION FORMAT 
DELTAX step size - x direction BortO. 3 
EPSI exponent for body radius 

EPSIUFe tolerance on velocity for a streamline call 

EPSIC=€ | tolerance on ali species for a streamline call 

EPSIH=€,, tolerance on total enthalpy for streamline call 

EPSIT=€,, tolerance on temperature iteration 

FH (N) total enthalpy profile, N = 1,2...LP2 

FK factor in ve 
*CPE free stream specific heat at constant pressure BetOsS 
FL characteristic length to dimensionalize 

species production term 

FSHE free stream static enthalpy 

GAMM ratio of specific heat, Cp/cr 

PE free stream pressure 

RHOE free stream density 

RHOSTG stagnation density 

CE (i) free stream species, I = 1,2...7 

Ry nose radius 

UINF free stream velocity 

ULOLIM low limit for velocity in streamline 

XS x starting value 
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SYMBOL DEFINITION FORMAT 
Z1s €, coordinate starting value E 10.3 
22S C2 coordinate starting value 
Z1L €, coordinate entire body 
u(N} velocity profile, N = 1,2,...LP2 
*FA coefficient for FKPSI in pre E 10.3 
FB FKPSI = K(N) (see p.8, Section ii) 
FC 
TESTRA test ratio for establishing turbulent viscosity 
model 
AB, law of the wall - E ve 


dimensional coordinate y 


ABs law of the wake - 
dimensional coordinate, Ye 
T : ; ; 
ABs laminar viscosity 
T : : : ; 
AB, eddy viscosity fit Indicators and 
dimensional coordinate y coefficients for 
T invoking turbulent 
ABg law of the wall - viscosity model 
streamline coordinate, ¥ components 
T (see p.7, Section 4,b) 
ABs law of the wake - 
streamline coordinate, VY, 
T ; : ; 
ABy eddy viscosity fit 
streamline coordinate, W 
T ; 
BR, law of the wall - 


dimensional coordinate, y 


PAS TTEOBUTTELITTTISTTD TT PET STEEP RET ET TSIPEULS PEs pester reels pee loeese HPS Mo bert ee er en La) TER ee TTF eI CET 


. As 


SYMBOL 


FLELIN 
FLETIN 
PRALIN 
PRATIN 
SCHLIN 
SCHTIN 


¥(N) 


*INDPRI 
NSR 
JINPUT 


TriPER 


DEFINITION 


law of the wake - 
dimensional coordinate, ve 


laminar viscosity 


eddy viscosity fit 
dimensional coordinate, y 


law of the wall - 
streamline coordinate, WV 


law of the wake - 
streamline coordinate, Y 


eddy viscosity fit 
streamline coordinate, W 


Lewis laminar number 
Lewis turbulent number 
Prandtl laminar number 
Prandtl turbuient number 
Schmidt laminar number 


Schmidt turoulent number 
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FORMAT 


E 10.3 


Indicators and 
coefficients for 
invoking turbulent 
viscosity model 
components 

(see p.7,Section 4,b) 


E 10.3 


physical coordinates for input profile 


normal to body 
print interval x direction 


# shock regions 


# times to bypass step size control halving 


% to be used to establish boundary layer je 6 aS) 


thickness 
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SYMBOL 
*ASR(N) 
BSR (N) 
CSR(N) / 
Z1R(N) 
Z2R(N) 
INDTYP (N) 
INDCOO(N) 


INDLAS (N) 


Note: ASR(N)-INDLAS(N) are repeated for each shock region, 


DEFINITION 


coefficients for th shock region geometry, 
N= 1, NSR 


tw pewndary. for yth shock region, N=1, NSR 
th i 

C2 boundary for N shock region, N =l, NSR 

geometric type, N=1, NSR 

end coordinate indicator, N=1l, NSR 


last shock region indicator, N=1, NSR 


N =1, NSR 


*INDP 
INDR 
INDLR 
INDLSR 
INDS5 
INDS6 
INDS7 
OGIVEH 
CGIVEK 
P(J) 


ac (K,I) 


A osarsnsen2n gaan sazasasanadisnonndiaantasinonsnenianniadaaiaen aliiiian ccasiadtdiiat tite tabasasiemeeeniiiaiia aati k et ek ell 7 


body subregion pressure type 
body subregicn geometric type 
last region indicator 

last subregion indicator 
spare 

spare 


spare 


coordinates of center of body arcs 
or ogives 


coefficients for pressure formulas, J =1,8 


coefficients for species mass fractions 
at wall, K = 1,5; I = 1,NS-1 


E 10.3 


he EFT 
eecetsenedttt oh 


Muaneaneriuenity 


HineoreEy 
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SYMBOL DEFINITION FORMAT 
HH (J) coefficients for total enthalpy at wall, J=1,5 

Ry body radius 

x, cumulative x limit for each subregion 

CONEAN cone angle 


Note: INDP-CONEAN are repeated for each body subregion 


* 
Indicates that this field is the first one on the input card, 
and subsequent fields until the next * are continuous through 
column 70 for as many cards as are necessary. 
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V. PROGRAM FEATURES 
A. Step Size Control 
The purpose of step size control is to allow the program 
to run as fast as possible within stability and truncation 
requirements. 
l. W direction 

Whenever the number of W intervals becomes twice 
the original input amount (as a result of swallowing) the program 
halves the number of points. This is of significance since the 
permissible stability control step size in the x direction varies 
with (AW?*. 

2. %X* direction 

The program takes two single steps and compares 
the velocity, species and enthalpy with the results of a single 
step of double size. Depending on input limit tolerances, the 
step size is doubled, unchanged or halved. However, the step 
size is always subjected to the stability control for maximum 
size. 

If a negative velocity, specie, enthalpy or temper- 
ature has been calculated, or depending on a sense switch being 
on, if a 5% out of monotonicity in u or H develops, the step size 
is halved up to seven times before the program stops. Provision 


has been made for overriding the halving for JINPUT times, so as 
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to allow zero species to fill in as a result of production and 
of diffusion in the WY direction. The above procedure for 
doubling, continuing and halving the step size which applies to 
the boundary layer also applies individually to u and C3 in the 
streamline calculations. 
B. Stability Control 

The expansion of the step size in the x direction which 
is described under "Step Size Control" is subsequently scrutinized 
to insure stability. 


4x is compared to Sf and accordingly as it is less 
n 


than (or equal to) or greater than, it is accepted or halved. 


iti rari 
Be LL n+1/2 n-1/2 
on Ee oe Ar B + Ay 2 


~ 


L —T aa 
eon : aL ae A al Ma+l/2 Ma-1/2 1 
mmnn}] * pu be 2 

xr 


ee = ny 
1g 
n n 
L L —T —T 7 
+ 
oO fF rrpiuiijaA one ok + Be preach tape “k 
= 3 3 
Cc mmnn gh! pu 2 st pe 
(e n Cc n 
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AY Av? AY? 
eon oo iin vale 


Ax must be S the minimum of in order 


to insure stability. 
Cc. Conversion from Physical to Streamline Coordinates 

1. The various input profiles are given as functions 
of the physical coordinate Y, which may be unevenly 
spaced. The number of W intervals equals the input 
number of Y intervals and is equal to L- l. 

2. The streamline coordinate WV = He pur‘ dy is calcu- 
lated and then equispaced over the thickness. 


3. The corresponding Y's are developed from 


a 
y=f z av with ¥(0) = 0. 
o pur 


4. A table lookup routine is used to linearly inter~- 
polate the dependent variables, u, Ci, H. 

5. The calculated value of ten is compared with the 
input value of 5 ae (see Input, Section IV.D). If 
the difference is greater than .1%, the following 
on line message is printed: "YCALC(L) = NACE) coe 

D. Continuation Procedure 
At any time during the running of the program sense 
switch 4 can be depressed and an on-line punch of the continua- 
tion quantities will occur. This makes it possible to continue 


the program at a later date. 


When the program is continued with on-line punched cards, 
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RESTAR must be set equal to zero. 


E. Output Control 
Based on the input quantity INDPRI, output at every 


a9) : 
n® station is obtained and based on NPSI, output at every m 


W level. There is also provision in the streamline calculation 


‘eR TT SRR TTT A nA TTT 


to obtaiil output whenever the step size doubles. 


F. Streamline Control 
An initial value or a complete body problem can be 


solved depending on whether we input a 0O or 1 for INDSTR. Thus 


there are two possible starting points for a streamline calcula- 


Pee rere 
m7 ETE TETPERT TT TEP PETTY PRET eTPETT PETES eT HTETTT TTT FeysTerrereeeerereysevrereee 


tion: 

| 1. The starting profile at the starting point on the 

| body. 

2. The starting profile at the point of intersection 

of the streamtube and the shock. 

| In addition, the first option can be additionally modi- 
fied tc constant edge conditions, i.e., Ug Et3) ae 


’ at = ’ s i es B 
Hey 3) He by setting INDSTR 1 


‘ at = F 
Cee 3) Cant 


“a 
G. 


Sense Switch Control 


Is available for diagnostic and monitoring purposes 


can be specifically determined by consulting the source language. 
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H. Additional Programming Logic : 
1. Momentum equation 
a. In order to step forward the required p 


is obtained as follows: 


n(w = H(y - 2 


a 


Using as our first approximation for T 


+ 
eae pet 0 a ~ Sy oa Biren = a Po mae 
fe) N NO not 
aa Py Paar : 
NO 2 NO p 
calculate 


ed 


C. C D D 
aE) Deg Gs sea ae / eyes) nee oe 
Joe kM 0 N 
+ 
Cuo | Px2tPo2 Snot DiotPoo 
2 = NOs oh Ii0 2 Dio 
™NO Not 


j = Og, Nz ,NO,NO* 


k=0,N,e 


Pe < 
Piet 
cycle at A. 
135 
foie 
rou 


we use T otherwise continue the iteration 


itl’ 
p(x) is chosen according to the indicator 


for the subregion and we compute 


Because the nature of the u profile is not 

*inear in the neighborhood of the wall, me is 
obtained by assuming u = av + bW and deter- 
mining a and b from ug and ug. The derivative at 
the 1/2 point is obtained by differentiating u. 

If a velocity is negative or if both sense switch 
5 is depressed and a 5% out of monotonicity occurs, 
then the step size is halved and the calculation 


repeated. If after halving seven times, the 


velocity is still negative, the program stops. 


2. Species equation 


If a species is negative, the step size is 


halved and the calculation repeated. If after 


halving the step size seven times, a species is still 


Lis uiasseball 
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negative, the program stops. oc is obtained by 
assuming C = Cy + a? + bY and determine a and b 
from Cz and Cs. The derivative at the 1/2 point 
is obtained by differentiating C. 
Enthaipy equation 

Halving ot the step size will result when 
a negative enthalpy occurs. 

If after seven halvings, this condition per- 
sists, the program will halt. 

Hy, is obtained by assuming H = H + av? + bY 
and determine a and b from Hz and Hg. The deriva- 


tive at the 1/2 point is obtained by differentiating 


He 


oaraseppsnnny teen 


= 
= 
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VI. PROGRAMMED STOPS AND PAUSES 


3 PROGRAM STOP NO, MEANING 
Displayed in 

Address Ficla of 

Storage Register 


MAIN Ly € <= 0 


77 Shock geometry is cone followed 
by parabola 


310 A. < 0 
1 
650 Kes 
s 
2575 LPs 2 99 
2655 L > 2 (LORIG) 
CED3E (J) 60 ir soeel 
71 Cron OG 1G aoe cin ernor 
; w,i 
RADBDY 550 € <0 
STEPSZ 603 (Pause) Ax has been halved 7 times 


Additional 7 times can be 
obtained by depressing START 


D 
UCSTRE 21 > «iS not in any shock region 
UEDGE (J) 30 ac oat 


VII. OPRRATING PROCEDURE 
A. Original Run 
Standard 7094 Fortran operating procedure is employed 
with RESTAR = i.0E0. Progam language is FORTRAN I1. 
B. Continuation Run 
Remove subregion parameter cards that were read in 
during previous run. Insert continuation cards after problem 
input and ahead of remaining subregion cards. Set RESTAR = OEO. 
There is present the option of modifying the original 
input or the continuation input by repunching the appropriate 
fields or cards. It should be noted that the continuation cards 


are necessarily in octal in order to preserve full significance. 


Ba 


j 


SST TTT TPIT , cole 
eerste tear ret ern TTT ETT 


APPENDIX 1 


DIMENSIONS OF VARIABLES 


Toasters itil 


Fiadiapivecgas ) cede 


APPENDIX 1 


tain 


Ez 
Ez 


DIMENSZTONS OF VARIABLES 


The units of the variables which appear in the equations 
are listed below. The conversion constants required for con- 


sistency are indicated. 


A. Momentum Equation 


1. The units for the quantities involved are 


os lb _sec* a slugs 
ee Ete 
u = ft/sec 
p = lbs/ft® 
dp _ 


ae lbs/ft*, s is nondimensional 


ere A : 
ry lbs/fit™, x is in feet 


rT = lbs/ft® 
he lbs sec 

Lt fte 
On ay cee 


lbs_sec* 


te R should 


2. In srder for P to have the units 


be expressed in the following units: 


Secret Ucn eae Detar ne pako/ines 


ft* RT DC/M _. mole 
R-K =——— 


gm 


2 
= 
Re gm ft 


[e) 
K mole sec” 


-i 


8.31657x10° joule {kilogram mole *) °K 


i 


Se eee tees gm 2 
ST As6: souleS 453.6 4 32.2 ft/sec 


m £t* 


7 
8.95805 x10 K (kilogram mole) sec 


m ft” 
K (gram mole) sec 


8.95805 x 10* 


3. In order for T to have the units °K the following 


adjustments must be made: 


D 
ees ft* [=] 4 
m ° 


T ™o sec art ~ ev sec* °K 
co _—_—_——— Cos = <i 
C £e part gm ft 
P sec K 


ait 
i 
LES) 


sesegan eae ott4e 


Taking the second term 


presronavenennye 


ev_sec* en eeu aes: 1 -ftllb 1 joules 
: pact om £t* *° lb ~“** sec” 1.356 joule 6.24x 10! av 
3 
= 1726x110 °° (°x) 
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Thus for all terms, this factor should be used, i.e. 
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D ore D 
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B. Species Equation 


1. The units for the quantities involved are 


= lb sec 
f£t* 


Se 


which is derived from 


se 
| 

xe) 

oc 


2se 4D should be of the same units as p for Eq. (7a), 


DO 
p. 5 of GASL TR-246. Hence, 
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Pro ~ article a SEE =: = x) 
amr ee (ft lbs sec)~ P 


&/, 


a gm ieee 
(ft lbs)¥2 sec*part2 | 453.6 gm | 


~11 lbs sec” -- 
ft* DO 


3.87871 x10 


Similarly for Pon! P oxo! Por 
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3. 0, ¥, K are dimensionless 


4. <A, B, C, K develop as follows (example for A shown below) 


£t fea Sai 1b sec® Lb sec’ 
gm/mole |_gm_ part*sec Bee: ft* 
mole 


2 (ee 
BiB aaa 
eee glee Few de [453.6 <a] AB bat i a] [é-0251 x 10** ok 
Et‘ gin* part lb 2.54 cm 12 in mole 


mole sec — 


= 4.9778 x10"? 
part 


5. D, E, F develop as follows (example for D shown below) 


ft em® lb sec® 
gm_ part sec ft* 


mole 


3 3 
D= mole cm™_sec lb 453.6 2 1 ES} 
rt? gm part lb {30.48 cm 


-2 mole sec — 


= 1.6018 x 10 pareicle 


6. Equation (9), p. 7 of GASL TR-246 would yield (neglect- 
M 
ing —°) 
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u Wa aa but “Fre are required hence 
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Sane 4 32.2 ft/sec” 6.0251 nase 16 = ft 
ne 25 2 Ib sec = 
W 1.94008 x 10 2 | Fee We 
Poo 
in Ay Bi, etc. yr should be dimensionless. This can 
be accomplished by setting 
ev 
D a : 
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! D \ 
a 19 Q2 
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Energy Equation 
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The units for the quantities involved are: 
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